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Long-term variations of solar activity, including the Grand minima, are 
believed to result from temporal variations of dynamo parameters. The simplest 
approximation of dynamo waves is applied to show that cyclic variations 
of the parameters can lead to an exponential growth or decay of magnetic 
oscillations depending on the variations frequency. There is no parametric 
resonance in a dynamo, however: the selective sensitivity to distinct frequencies, 
characteristic of resonant phenomena, is absent. A qualitative explanation for this 
Ending is suggested. Nonlinear analysis of dynamo-waves reveals the hysteresis 
phenomenon found earlier in more advanced models. However, the simplified 
model allows a computation of a sufficiently large number of dynamo-cycles 
for constructing the distribution function of their amplitudes to reproduce 
qualitatively two modes of solar activity inferred recently from cosmogenic 
isotope content in natural archives. 
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INTRODUCTION 


Hydrodynamic flows generating cosmic magnetic fields usually vary with time. Parame¬ 
ters of the generation processes - hydromagnetic dynamos - are, therefore, time-dependent. 
Variations of dynamo-parameters can be either periodic, as in the case of modulation by 
density waves in galactic dynamos, or irregular, as in stellar convection zones. The case 
of irregular variations currently attracts much attention in relation to the Grand minima 
of solar activity (cf., e.g., the review of Charbonneau 2010), events similar to the famous 
Maunder minimum when there were almost no sunspots for about 70 years (Hoyt & 
Schatten 1996). 

The number of large-scale convective cells generating magnetic fields in the Sun is 
not large (Gilman & Miller 1985) so the longitude-averaged dynamo parameters fluctuate 
considerably (Olemskoy et al. 2013). Fluctuations of dynamo-parameters are believed to 
be responsible for the observed variations of solar cycle durations and amplitudes and even 
for the events of Grand minima. Dynamo models with fluctuating parameters generally 
reproduce these observed phenomena (cf., e.g.. Moss et al. 2008; Usoskin et al. 2009; 
Ghoudhuri & Karak 2012; Karak & Ghoudhuri 2013; Olemskoy & Kitchatinov 2013). 
Relations between the fluctuation properties and the resulting variations in magnetic 
cycles remain, however, unknown. In particular, the important theoretical question on 
the possibility of parametric resonance in dynamoe 3 remains unanswered. It is not clear 
whether there are special frequencies in the fluctuations spectrum to which the field 
generation processes are especially sensitive. In some dynamo models, the characteristic 
fluctuation time was put equal to a half of the magnetic cycle period, which corresponds to 
the parametric resonance condition for harmonic oscillator, but an enhanced reaction on 
fluctuations of such a duration was not found. The reason probably is that contemporary 
dynamo models are numerical and it is difficult to find functional relations from the results 
of numerical computations. 

Probably for the same reason, the publications on resonant phenomena in dynamos are 
not many. The dynamo equations do not possess sources of magnetic fields. The resonance 
problem, therefore, does not have an analogue of driving force, and only parametric 

tThe non-periodic character of the parametric fluctuations is significant but not of principal 
importance. Stochastic parametric resonance (Klyatskin 1980) should be considered for this case. 
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resonance can be considered. Schmitt & Rudiger (1992) and Moss (1996) studied the 
parametric resonance for galactic dynamos. Gilman & Dikpati (2011) and Moss & Sokoloff 
(2013) considered the resonance phenomena for stellar convective dynamos. An interesting 
version of parametric resonance in nonlinear dynamos was found by Reshetnyak (2010); 
the resonance can result from a time-delay in the reaction of dynamo-parameters to the 
generated magnetic held. 

This paper suggests using a simplihed dynamo model to treat the resonance problem 
analytically. Analytical solutions may be helpful in understanding the results of more 
complicated numerical models. Such an approach proved to be productive in studies 
of magneto-rotational instability where a simplihed analytical model revealed similarity 
rules conhrmed subsequently by numerical computations (Kitchatinov & Rudiger 2004). 
Sufficient simplihcation for an analytical treatment of the resonance problem give an 
(local) approximation of dynamo-waves. In this approximation, the dynamo-wave length 
is assumed short compared to the scale of variations of dynamo parameters. Such a scale 
ratio does not hold for the Sun and for corresponding numerical models of global dynamos. 
However, it is difficult to hnd a numerical model whose results were not interpretable in 
terms of the Yoshimura (1975) law for dynamo-waves propagation. 

We shall see that the parametric modulation of dynamo waves can lead to an exponential 
growth or decay of the wave amplitude depending on the frequency 7 of variation of the a- 
parameter of magnetic held generation. This modulation, however, is not a resonance. The 
enhanced selective reaction to special frequencies, characteristic of resonant phenomena, is 
not present in dynamos. A qualitative explanation for this circumstance will be suggested. 

Another problem where the dynamo-wave approximation can be useful is related to 
the two distinct modes of solar activity inferred recently by Usoskin et ah (2014) from 
radiocarbon data. The two modes may be related to the hysteresis in nonlinear dynamos: 
allowance for the dependence of turbulent dihusivity on the magnetic held reveals the 
coexistence of two stable dynamo solutions with diherent amplitudes of magnetic cycles 
that are realised depending on the pre-history of dynamo parameter variations (Kitchatinov 
& Olemskoy 2010). Dynamo-hysteresis was also found in 3D numerical simulations of 
Karak et al. (2015). However, a moderate number of magnetic cycles covered by computations 
with complicated models hinders their statistical analysis. We shall see that allowance for 
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nonlinearities in the dynamo-wave approximation also shows the hysteresis phenomenon. 
Millions of magnetic cycles can be computed with this approximation. The distribution 
function of the amplitudes of these cycles has two maxima similar to the dual distribution 
found by Usoskin et ah (2014) from historical data on solar activity. 


BASIC EQUATIONS 


We start from the mean-field dynamo equation 


dB 


rot {aB + V x B — rj rot B) 


( 1 ) 


(cf., e.g., Zeldovich et ah, 1983). The velocity V corresponds to the shear flow of, e.g., 
non-uniform rotation. Cartesian coordinates are used with the y-axis aligned with the 
velocity vector V and the x-axis along the velocity gradient. 


V = (0, Sx, 0). 


( 2 ) 


Coefficients in the equation ([T]) and the shear value S are assumed position-independent, 
i.e. the scale of magnetic field variation is assumed small compared to the characteristic 
scale of the fluid inhomogeneity (short dynamo-waves). For the convection zone of a star, 
the |/-axis of our coordinate system corresponds to the azimuthal direction, the x-axis 
points in the direction of the angular velocity gradient and the ?/ 2 ;-plane is the isorotation 
surface. The magnetic field, 

B=(-|4.B.o). (3) 

includes toroidal (B) and poloidal {—dA/dz) components. 

The field B and the poloidal field potential A depend only on the coordinate 2 ; and 
time. The equation system for the two field components reads 


dB 

dA 

dt 


dA d^B 


aB 


d^A 


( 4 ) 


where the contribution of the a-effect into the toroidal field generation is omitted (the 
ahl-dynamo). 

We consider wave-type solutions of A, i? oc exp{ikz) and normalize the equations (jlj) 
to dimensionless units. Time is measured in units of (fc^r^o)”^. The t]q is the characteristic 
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value of diffusivity, rj = rjQfj {t) ~ 1 ), and similar scaling is used for the a-coefficient: 
a = aoa. The potential A is measured in units of aoBQ/{r]Qk‘^)] Bq is the characteristic 
value of the held B and simultaneously the unit for held measurement. We keep the same 
notations for normalised variables as used before for their dimensional counterparts to 
hnd 


B = —iVA — fjB, 

A = aB — rjA, (5) 


where 


V = 


aoS 


( 6 ) 


is the dynamo-number. The parameters a. and fj may depend on time explicitly, as they 
do in the linear problem on parametric resonance, or implicitly via their dependence on 
the magnetic held amplitude in a nonlinear problem. 


ABSENCE OF PARAMETRIC RESONANCE 

If the parameters of Eqs. ([5]) do not depend on time (17 = d = 1), the threshold value of 
the dynamo-number for emergence of non-decaying oscillations V = 2 (the amplitude of 
the oscillations grows exponentially with time for the larger dynamo number and decays 
with time for the smaller D). The eigenvalue problem {A,B oc exp(—iwt)) gives two 
solutions for this case, 

ui = 1, U2 = -I - 2i, (7) 

the hrst of which corresponds to non-decaying oscillations. Let us consider how this 
solution changes with allowance for small variations of the d-parameter with time: 

d = 1 -I- ea{t), ( 8 ) 

where e -C 1 is the amplitude of the parameter perturbations and a ~ 1. On eliminating 
A from equations (j5]), we hnd 

B + 2B + {l + 2i + 2iea{t))B = 0. (9) 

The problem of parametric resonance concerns periodic variations of the parameter 
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with time: a = sin( 7 t). We search for the solution in the form 

OO 

S = e‘"‘ ^ (10) 

n=—oo 

It can be seen that a decomposition of the coefficients bn in powers of the small parameter 
£ starts with With accuracy up to the second order in e, we find the amplitudes of 
‘satelite’ oscillations, 


bi 


g bp 

7 ( 7 - 2 ^- 2 )’ 


g bp 

7(7 + 2 i + 2 ) ’ 


and the eigenvalue 


a = s + iAu, 


e '^ 7^-8 
¥74 + 64’ 


Au 


7 ^ + 8 
¥74 + 64’ 


( 11 ) 


( 12 ) 


including the (real) growth rate s and frequency shift Au. The approximate solution 
(1TT|) . flT^ is valid for the small ratio e/'j 1 . It can be seen from the equation (IT^ for 
s that either exponential growth (for 7 > \/8) or decay (for 7 < -s/S) of the oscillations is 
possible. 

It is of course possible to solve Eq. fllUp numerically to confront the computations 
with our approximate analytical solution. Figure 1 gives examples of such a comparison 
for decaying and growing oscillations. The initial conditions for the computations were 
B = 1 and B = —i (for f = 0). Equations flTT]) and flT^ show that the main consequence of 
the parametric modulation - of the first order in e - is the deformation of the oscillations’ 
shape. Variation of the oscillation amplitude is the second-order effect. Part (a) of the 
Fig. 1 clearly shows the oscillation asymmetry: negative deflections have larger amplitude. 

The oscillations grow exponentially for sufficiently large frequency 7 . The growth does 
not mean a resonance, however. The dependence of the growth rate s on the modulation 
frequency 7 of Eq. flT^ is smooth and does not show any selected frequencies characteristic 
of the resonant phenomena. There is, therefore, no parametric resonance in the dynamo. 
The reason for this can be seen by comparison of the equation ([9]) with the Mathieu 
equation for the parametric resonance, which for further discussion is convenient to write 
as 

X + oj^X = £sin( 7 f)X. (13) 


With parameter oscillations omitted {e = 0), there are two eigenfrequencies u = ituo. 
Following the simplest perturbation method (not practical but convenient for qualitative 
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Fig. 1. Damped oscillations for 7 = 1, e = 0.3 (a) and growing oscillations for 7 = 5, e = 0.7 
(b). The full line shows the numerical solution. The dotted line shows the analytic trends 
±exp(st) for the oscillation amplitude corresponding to the Eq. (11211 . 


analysis), we can search for the first-order corrections in e: by snbstitnting the ‘nnpertnrbed’ 
solntion X oc exp(±za;ot) into the right-hand side of (IT^ . This gives a ‘driving force’ 
oscillating with the freqnencies icuo ± 7 (part (a) in Fig. 2). For 7 = 2a;o, the freqnencies 
—uq + 7 and cjo — 7 nieet the eigenfreqnencies of the nnpertnrbed system to prodnce a 
resonance. This is why the growth rate for parametric resonance is of the hrst order in e 
(Landan & Lifshitz 1969). The case with dynamo is different (part (b) in Fig. 2). There 
is only one eigenvalne (jT]) for non-damped oscillations. Displacements to the left or to the 
right by the valne 7 on the freqnency axis from this eigenfreqnency, shown by the npper 
bows in the part (b) of the Fignre, do not lead to other eigenfreqnencies. The next step of 
the pertnrbation techniqne shonld be nndertaken to get a hnite growth rate. This second 
step gives once again displacements to the left and right on the freqnency axis shown by 
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(a) 

— U<’() 0 a,'(| 

(b) 


Fig. 2. (a) Schematic representation of the origin of parametric resonance for oscillations 
governed by the Mathieu equation (b) Non-resonant parametric modulation in a 

dynamo (cf. explanation in the text). 

the lower bows in Fig. 2(b). The second step of the perturbation technique returns to the 
original eigenfrequency to give a hnite growth rate of the order fll 2 p . This happens 
for any modulation frequency 7 . Parametric modulation in a dynamo is, therefore, ‘not 
selective’ in 7 , i.e. it is not resonant. Parametric resonance is probably possible only for 
systems with not less than two eigenfrequencies of non-damped oscillations, and a resonant 
frequency equal to the difference between these eigenfrequencies. 

We turn now to the case of a{t) in equation ([9j) being a random function of time. 
For this case, the equation has to be solved numerically. Random function with known 
statistical properties can be realized with the equations 

da a ai dai 

dt r r ’ dt 

where r is the correlation time and 5 ^ is a given random forcing. In the computations to 
follow, g = g^JdrjAt, where At is the numerical time step and g is the random number 
with normal (Gaussian) distribution, the mean value equal to zero and RMS value equal 
to one. The number is renewed on each time step independent of its preceding value. The 
function g represents a short-correlated random process but a{t) is a smooth function of 
time with zero mean. For short At -C r, it represents a stationary (in the statistical sense) 
random process with the correlation function (Olemskoy & Kitchatinov 2013) 

{a{to -b t)a{to)) = {\t\ /t + 1) exp(- | t | /r), (15) 

where the angular brackets signify averaging over time to (see also Rempel 2005). 

Figure 3 shows a characteristic example of the oscillations resulting from a joined 
solution of Eqs. (E]) and 0141) . Random wandering of the oscillation amplitude away from 


ai g 

T T 


(14) 


jj’ 1 ~ 
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Fig. 3. An example of linear oscillations in the presence of random fluctuations of the a-effect. 

It results from solving jointly the equations ([9]) and (|14p for e = 0.3 and r = 1. 

its initial value is typical of linear oscillation with stochastically varying parameters 
(Klyatskin 1980). Figure 3 also shows this type of behaviour. The diffusive wandering of 
the dynamo-wave amplitude is by all probability not realistic. It is stabilized by nonlinear 
effects, to the account of which we now turn. 

HYSTERESIS IN NONLINEAR DYNAMOS 

Dependence of the coefficients of the equation ([5]) on the magnetic field is now allowed 
for. The dynamo-wave approximation can still be applied to the nonlinear problem if the 
coefficients a and fj are the functions of the amplitude f3 =\ B \ of the wave (see Noyes et 
ah 1984 for more details). The equation system ([5]) is no longer reducible to one equation 
similar to Eq. ([9]) in the nonlinear problem. 

Instead of the Eq. ([H]), we now have for the d-coefficient 

4 = (1 + ^ (l - - h^arctg(«) . (16) 

where the dependence on the magnetic field is prescribed after the prediction of the quasi- 
linear theory of the a-effect (Rudiger & Kitchatinov 1993). Magnetic quenching of the 
a-effect is usually the only nonlinearity in dynamo models. It suffices for stabilizing the 
filed growth. The d value decreases with increasing f3. The effective dynamo number, 
considered as a function of the magnetic field, reduces to its critical value to stabilize the 
field growth. 
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However, not only the a-effect but the eddy diffusion as well can depend on the 


magnetic field. Quasi-linear theory of turbulent transport coefficients predicts 


V = 


8/32 


1 + 


4 + 


(1 + /? 


2\2 


+ 


/32 -5 


arctg(/9) 


(17) 


for the dependence (Kitchatinov et ah 1994). The fj, as well as a, decreases steadily with (3. 
The effective dynamo number Pgff = 'Da/rf', however, is no longer a monotonic function 
of the magnetic field. For large /3 S> 1, oc and the nonlinear saturation of the 
field growth still takes place. For small /9 -C 1, however, the effective dynamo number 
increases with the field strength. This leads to the hysteresis phenomenon. 



3 ) 

Fig. 4. Illustration of hysteresis in nonlinear dynamos. In the range of 0.6 < T> < 2, the 
equations ([5]) with the coefficients (|16p and (|17l) have two stable solutions of i) oscillations 
with a steady and relatively large amplitude j3 > 1, and ii) decaying oscillations of small 
amplitude. The dotted line shows the watershed between the ‘attraction regions’ of these 
two solutions. 


In the absence of a-effect fluctuations (e: = 0 in Eq. (lT6|) h there are two stable solutions 
of Eqs. (jS]) in the range of the dynamo numbers 0.6 < T> < 2 realised depending on the 
initial condition. Figure 4 illustrates this situation. Computations for a series of increasing 
dynamo numbers and a small initial amplitude [3 show decaying fields up to dynamo 
number T> = 2. When the number V reaches its critical value oi T>c = 2, the decay is 
replaced by field growth. This leads to an increase in the effective dynamo-number and the 
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amplitude of the eventually steady oscillations jumps to a finite value of /? > 1. A further 
increase in dynamo number results in a smooth increase in the amplitude. If the dynamo 
number is then changed to a decrease and the amplitude of the preceding computation 
is taken as the initial condition, oscillations of finite amplitude survive for P < 2 up to 
V = 0.6 but decay indefinitely for still smaller T> < 0.6. The equation Defr(/d) = 2 for 
the amplitude of the dynamo-wave has two solutions in the range of 0.6 < V < 2. The 
solution of the relatively small amplitude is however unstable. This solution is shown by 
the dotted line in Fig. 4. This line is the boundary between the ‘attraction regions’ of the 
oscillations with steady and relatively large amplitude (full line) and decaying oscillations. 
Dynamo waves decay if their initial amplitude is below the dotted line. Above this line, 
the amplitude grows initially and then stabilizes. 

Fluctuations in dynamo parameters can provoke irregular transitions between the two 
solutions. This possibility attracts attention in connection with the two modes of solar 
activity inferred by Usoskin et ah (2014) from the content of cosmogenic nuclides in natural 
archives. The distribution function of the amplitudes of solar cycles for about 11 thousand 
years that they obtained has two maxima. Such a dual distribution was interpreted as a 
manifestation of two distinct modes of solar activity. 

Recently Karak et ah (2015) found the dynamo-hysteresis in 3D numerical experiments. 
Their computations near the critical dynamo number show an irregular intermittency of 
magnetic cycles of relatively large amplitude with weak-field epochs. However, the time 
consuming 3D computations did not provide sufficiently vast statistics of magnetic cycles 
for finding their distribution function. This can be done with the simplified dynamo- 
equations of this paper. Figure 5 shows the distribution function for one million magnetic 
cycles computed by solving numerically equation ([5]) with allowance for the nonlinearities 
of equations (ITB]) and flTT)) and the a-effect fluctuations of model flTT]) . The dual distribution 
function is qualitatively similar to that found by Usoskin et ah (2014) from historical data 
on solar activity. 

It may be noted that the dynamo-hysteresis and irregular transitions between dynamo 
modes with different amplitudes of magnetic cycles are not restricted to the special case of 
the nonlinearities of Eqs. ([T6|) and fITTj) . The hysteresis is present whenever the dependence 
of the ratio a/if on the magnetic field has a maximum, e.g., for the simplest model of 
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1 10 
Cycle amplitude 


Fig. 5. Distribution function for the amplitudes of 10® cycles obtained by solving numerically 
equation JS]) with allowance for the nonlinearities of Eqs. (|lb|) . (I17p and fluctuations of the 
a-effect. Computations were performed for P = 2.1, e = 1 and r = 1. 

d = (1 + and 17 = (1 + / 5 )“^. 


CONCLUSION 

The simplest approximation of the short dynamo-waves of this paper certainly does not 
allow quantitative reproduction of observations of solar activity. This should be the subject 
of detailed numerical models. However, the approximation allows analytical solutions, 
which can be useful in interpreting the results of the numerical computations. Variations 
of solar activity, including such impressive phenomena as its Grand minima, are related 
to fluctuations in dynamo parameters. The majority of models for the Grand minima 
assume the characteristic time of the fluctuations equal to the duration of the (11-year) 
activity cycle, i.e. to a half of the entire 22-year magnetic cycle that corresponds to the 
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doubled frequency of the magnetic cycle. An enhanced response to fluctuations of such 
duration as expected for parametric resonance was not found however. Solution of the 
parametric modulation problem for dynamo-waves shows that parametric resonance is 
not present in dynamos. Exponential growth or decay of the wave amplitude depending 
on the frequency of the parametric variations is possible but a selective sensitivity to 
special frequencies characteristic of the resonant phenomena does not come out. This is 
related to peculiarities of the dynamo mechanism (cf. Fig. 2 and the discussion around). 
It becomes clear why the Grand minima statistics can be reproduced with either durable 
(Moss et al. 2008) or rapidly varying (Olemskoy et ah 2013) fluctuations. The fluctuation 
parameters should be inferred from observations. 

The dynamo-wave approximation turns out to be useful also in the problem of the 
origin of the two modes of solar activity indicated by the new data (Usoskin et al. 2014). 
The two modes may be explained by the hysteresis in nonlinear dynamos ( Kitchatinov 
& Olemskoy 2010; Karak et al. 2015). Hysteresis is also found in the dynamo-wave 
approximation, which permits its detailed analysis (Fig. 4). In contrast to more complicated 
models, computations with this approximation are not numerically demanding. It is 
therefore possible to accumulate the statistics of computed magnetic cycles sufficient 
for constructing a distribution function of the cycles’ amplitudes. A dual distribution is 
obtained qualitatively similar to that found by Usoskin et al. (2014) from radiocarbon 
data on solar activity. 

This work was supported by the Russian Foundation for Basic Research (Project M- 
13-02-00277). 
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